MODULE qcim_geometry

   USE qcim_global_consts
   USE qcim_global_types

   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: qcim_define_geometry

!FIXME dirty hack to allow for plate potentials in geom.qcim input file
REAL(REALK), SAVE :: plate_pot
CHARACTER(1), SAVE :: field_or_potential  ! 1==Field, 2==Potential

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE qcim_get_geometry(fem_order,ngrid,origin,lambda, &
                                use_images,pshift)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: fem_order
      INTEGER(INTK), INTENT(OUT) :: ngrid(3)
      REAL(REALK),   INTENT(OUT) :: origin(3), lambda(3)
      INTEGER(INTK), INTENT(OUT) :: use_images
      REAL(REALK),   INTENT(OUT) :: pshift

      REAL(REALK)  :: hi(3), field_input
      CHARACTER(9) :: use_images_char

      OPEN(UNIT=IO2, FILE='geom.qcim', STATUS='OLD',  &
           ACCESS='SEQUENTIAL', FORM='FORMATTED')
      REWIND(IO2)

      READ(IO2,'(F20.15)') origin(1)
      READ(IO2,'(F20.15)') origin(2)
      READ(IO2,'(F20.15)') origin(3)
      READ(IO2,'(F20.15)') hi(1)
      READ(IO2,'(F20.15)') hi(2)
      READ(IO2,'(F20.15)') hi(3)
      READ(IO2,'(F20.15)') lambda(1)
!FIXME allow for non-cubic FEs
      lambda(2:3) = lambda(1)

      ! Set potential on z-plates
      ! Note potential at z=0 is arranged to be 0
      READ(IO2,'(A)') field_or_potential
      IF (field_or_potential == 'F') THEN
         READ(IO2,'(F20.15)') field_input
         plate_pot = - hi(3) * field_input
      ELSE IF (field_or_potential == 'P') THEN
         READ(IO2,'(F20.15)') plate_pot
      ELSE 
         STOP 'ERROR: Need to input field or potential on plate'
      END IF

      ! Flag for whether we have image charges in Laplace equation BC
      READ(IO2,'(A)') use_images_char
      IF (use_images_char == 'NO_IMAGES') THEN
         use_images = 0
      ELSE IF (use_images_char == 'USEIMAGES') THEN
         use_images = 1
      ELSE IF (use_images_char == 'GHOST') THEN
         use_images = 2
      ELSE
         STOP 'error in IMAGES input keyword'
      END IF
      ! Add constant shift to potential over all space
      READ(IO2,'(F20.15)') pshift

      ! Check that geometry is symmetric wrt z about z=0
      IF (hi(3) /= -origin(3) ) STOP 'ERROR: box is not centered about z=0'

      ! Number of elements must be integer, so input lambda has to be
      ! checked, and if necessary, adjusted

      ! Start with ideal ngrid
      ngrid(:)  = (hi(:) - origin(:)) / lambda(:) + 1
      ! Make multiple of fem_order + 1
      ngrid(:)  = ngrid(:) - MOD(ngrid(:), fem_order) + 1
      ! Adjust lambda to be consistent
      lambda(:) = (hi(:) - origin(:)) / (ngrid(:) -1)

      ! Since lambda must be same for all x,y,z, we have to adjust x,y lengths
      lambda(1:2) = lambda(3)
      ! Number of grid points for new x-y lambda
      ngrid(1:2)  = (hi(1:2) - origin(1:2)) / lambda(1:2) + 1
      ! Make multiple of fem_order + 1
      ngrid(1:2)  = ngrid(1:2) - MOD(ngrid(1:2), fem_order) + 1

!      ngrid(1:2)  = (hi(1:2) - origin(1:2)) / lambda(1:2)
!      ngrid(1:2)  = ngrid(1:2) - MOD(ngrid(1:2),fem_order) + 4

      CLOSE(UNIT=IO2, STATUS='KEEP')

      IF ( PRODUCT(ngrid(1:3)) == 0 ) STOP 'ERROR: grid-spacing too large'

   END SUBROUTINE qcim_get_geometry

!-------------------------------------------------------------------------------

   SUBROUTINE qcim_write_boundary(ngrid,origin,lambda)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: ngrid(3)
      REAL(REALK),   INTENT(IN) :: origin(3), lambda(3)

      INTEGER(INTK) :: ix, iy, iz, n
      REAL(REALK) :: tmp
      LOGICAL :: nb(3)

      OPEN(UNIT=IO2, FILE='bpot.qcim', STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)

      WRITE(IO2) PRODUCT(ngrid(1:3)) - PRODUCT(ngrid(1:3)-2)
      n = 0
      DO iz = 1, ngrid(3)
         nb(3) = (iz /= 1) .AND. (iz /= ngrid(3))
         DO iy = 1, ngrid(2)
            nb(2) = (iy /= 1) .AND. (iy /= ngrid(2))
            DO ix = 1, ngrid(1)
               nb(1) = (ix /= 1) .AND. (ix /= ngrid(1))
               IF (ALL(nb(:))) CYCLE
               ! Boundary points only
               tmp = qcim_get_bc_value(origin,lambda,ngrid,ix,iy,iz)
               WRITE(IO2) ix, iy, iz, tmp
               n = n + 1
            END DO
         END DO
      END DO
      IF (n /= (PRODUCT(ngrid(1:3)) - PRODUCT(ngrid(1:3)-2))) STOP 'ERROR: bpot.qcim'
      WRITE(IO1,*) "Writing fixed boundary potential to 'bpot.qcim'..."
      WRITE(IO1,*) "nboundary pts =", n
      
      CLOSE(UNIT=IO2, STATUS='KEEP')

   END SUBROUTINE qcim_write_boundary

!------------------------------------------------------------------------------

   FUNCTION qcim_get_bc_value(origin,lambda,ngrid,ix,iy,iz)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN)  :: origin(3), lambda(3)
      INTEGER(INTK), INTENT(IN)  :: ngrid(3)
      INTEGER(INTK), INTENT(IN)  :: ix, iy, iz
      REAL(REALK) qcim_get_bc_value

      REAL(REALK) :: tmp, xyz(3)
!      REAL(REALK), PARAMETER :: xyz1(3) = (/20d0,20d0,20d0/), Q1 = 5.0d0
!      REAL(REALK), PARAMETER :: xyz2(3) = (/15d0,15d0,15d0/), Q2 = 5.0d0
!      REAL(REALK), PARAMETER :: xyz3(3) = (/-20d0,-20d0,-20d0/), Q3 = 5.0d0

!      CALL RANDOM_NUMBER(tmp)

!---------
!      xyz(:) = origin(:) + lambda * (/ix-1,iy-1,iz-1/)
!      tmp = 0d0 
!      ! 1 charge
!      tmp = tmp + Q1/SQRT(DOT_PRODUCT(xyz-xyz1,xyz-xyz1))
!      ! 2 charge
!      tmp = tmp + Q2/SQRT(DOT_PRODUCT(xyz-xyz2,xyz-xyz2))
!      ! 3 charge
!      tmp = tmp + Q3/SQRT(DOT_PRODUCT(xyz-xyz3,xyz-xyz3))
!---------

!      tmp = 1d0

!---------

      ! Two XY plates
      IF (iz == 1) THEN
         ! xy-face, z-min
         tmp = -plate_pot
      ELSE IF (iz == ngrid(3)) THEN
         ! xy-face, z-max
         tmp = plate_pot
      ELSE
         ! 4 other faces use linear interpolation of plate potentials
         tmp = -plate_pot + 2d0 * plate_pot * ( (DBLE(iz)-1) / (DBLE(ngrid(3))-1) )
      END IF

!---------

      qcim_get_bc_value = tmp

   END FUNCTION qcim_get_bc_value

!------------------------------------------------------------------------------

   SUBROUTINE qcim_define_geometry(fem_order,ngrid,origin,lambda, &
                                   use_images,pshift)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: fem_order
      INTEGER(INTK), INTENT(OUT) :: ngrid(3)
      REAL(REALK),   INTENT(OUT) :: origin(3), lambda(3)
      INTEGER(INTK), INTENT(OUT) :: use_images
      REAL(REALK),   INTENT(OUT) :: pshift

      ! Interface boundary geometry
      CALL qcim_get_geometry(fem_order,ngrid,origin,lambda,use_images,pshift)

      ! Write QCIM boundary potential to file
      CALL qcim_write_boundary(ngrid,origin,lambda)

   END SUBROUTINE qcim_define_geometry

!------------------------------------------------------------------------------

END MODULE qcim_geometry

